df <- data.frame(year = 1994:1998,
born = c(69666, 69771, 67638, 67648, 66174))
ggplot(df, aes(x = year,
y = born)) +
geom_bar()
ggplot(df, aes(x = year,
y = born)) +
geom_bar(stat = "identity")
ggplot(df, aes(x = year,
y = born)) +
geom_bar(stat = "identity") +
ylab("Antal fødte")+
xlab("År")
11000/22
10500/22
library(MCMCpack)
install.packages("MCMCpack", dependencies = TRUE)
library(MCMCpack)
MCbinomialbeta(0.32*536, 536, alpha = 1, beta = 1, mc = 10000)
share <- MCbinomialbeta(0.32*536, 536, alpha = 1, beta = 1, mc = 10000)
summmary(share)
summary(share)
share1 <- MCbinomialbeta(0.32*536, 536, alpha = 1, beta = 1, mc = 10000)
share1 <- MCbinomialbeta(0.23*371, 371, alpha = 1, beta = 1, mc = 10000)
share1 <- MCbinomialbeta(0.32*536, 536, alpha = 1, beta = 1, mc = 10000)
share2 <- MCbinomialbeta(0.23*371, 371, alpha = 1, beta = 1, mc = 10000)
summary(cbind(share1, share2))
summary(share1)
summary(share2)
summary(share1 - share2)
prop.test(c(.32*536),536)
prop.test(c(.32*536, .23*371),c(536, 371))
prop.test(c(.32*371),371)
2 * 3 * 3 * 3
625*24
1 7 0.00417
1 / 0.00417
240 / 15
install.packages("swirl") # install the package library(swirl) # load the package install_course_github("kosukeimai", "qss-swirl") # install the co
library(swirl)
rm(lost = ls())
rm(list = ls())
swirl()
5 + 7
bye()
swirl()
install.packages("swirl") # install the package library(swirl) # load the package install_course_github("kosukeimai", "qss-swirl") # install the co
install.packages("swirl")
install_course_github("kosukeimai", "qss-swirl")
library(swirl)
install_course_github("kosukeimai", "qss-swirl")
swirl()
8 - 2
sqrt(10)
10^2
sqrt(9)
result <- 8 - 2
result
course <- "social science"
course <- "learning R"
world.pop
world.pop[4]
world.pop[c(1, 4)]
7
length(world.pop)
seq(1950, 2010, 10)
year <- seq(1950, 2010, 10)
names(world.pop) <- year
x[1] <- NA
290.35/13
22*13
x <- "x"
str(x)
x <- 1:2
str(x)
x <- 1:100
x <- 1:100/100
library(swirl)
swirl()
summary(afghan)
32.39
4
table(afghan$province, afghan$employed)
tapply(afghan$employed, afghan$province, mean)
info()
table(afghan$employed, afghan$province)
tapply(afghan$province, afghan$employed, mean)
prop.table(table(afghan$employed, afghan$province))
sum(is.na(afghan$violent.exp.taliban))
matrx(1:10, nrow = 2, ncol = 5, byrow = TRUE)
matrix(1:10, nrow = 2, ncol = 5, byrow = TRUE)
x <- matrix(1:10, nrow = 2, ncol = 5, byrow = TRUE)
apply(x, 1, mean)
z
z[[1]]
z[["z3"]]
155000/24
155000/72
155000 / 4
15500 / 12
social <- read.csv("c:/users/jod.dbp/dropbox/undervisning/cbs/pmsd/lectures/prediction/social.csv")
social <- read.csv("~/dropbox/undervisning/cbs/pmsd/lectures/prediction/social.csv")
social$age <- 2008 - social$yearofbirth
social <- subset(social, messages == "Control" | messages == "Neighbors")
mod1 <- lm(100*primary2008 ~ messages + age, data = social)
mod1
predicted <-
predict(mod1, newdata = data.frame(age      = rep(c(22:75), 2),
messages = rep(c("Control", "Neighbors"), each = 54)))
plot(22:75, predicted[1:54], type = "n", ylim = c(10, 50), ylab = "Turnout", xlab = "Age", main = "Neighbor treatment over age")
lines(22:75, predicted[1:54], col = "Blue", lwd = 2)
lines(22:75, predicted[55:108], col = "Red", lwd = 2)
legend(x = 22, y = 50, legend = c("Neighbors", "Control"), lty = c(1,1), col = c("Red", "Blue"), lwd = 2, bty ="n")
lines(c(50,50), predicted[c(83, 93)])
lines(c(50,60), predicted[c(93, 93)], lty = 2)
text(y = predicted[93] - 2, x = 45, paste(10*round(coef(mod1)[3], 3), "(10*age) ="))
lines(c(60,60), predicted[c(39, 29)])
lines(c(50,60), predicted[c(29, 29)], lty = 2)
text(y = predicted[39] - 2, x = 65, paste("=", 10*round(coef(mod1)[3], 3), "(10*age)"))
lines(c(30,30), predicted[c(9, 63)])
text(x = 36, y = mean(predicted[c(9, 63)]) + 3, paste("=", round(coef(mod1)[2], 2), "(Neighbors)"))
summary(lm(100*primary2008 ~ age*messages, data = social))
mod1 <- lm(100*primary2008 ~ age*messages, data = social)
predicted <-
predict(mod1, newdata = data.frame(age      = rep(c(22:75), 2),
messages = rep(c("Control", "Neighbors"), each = 54)))
plot(22:75, predicted[1:54], type = "n", ylim = c(10, 50), ylab = "Turnout", xlab = "Age", main = "Neighbor treatment over age")
lines(22:75, predicted[1:54], col = "Blue", lwd = 2)
lines(22:75, predicted[55:108], col = "Red", lwd = 2)
legend(x = 22, y = 50, legend = c("Neighbors", "Control"), lty = c(1,1), col = c("Red", "Blue"), lwd = 2, bty ="n")
lines(c(25,25), predicted[c(4, 58)])
lines(c(70,70), predicted[c(49, 103)])
coef(mod1)
round(coef(mod1)[3] + 25*coef(mod1)[4], 2)
mod1 <- lm(100*primary2008 ~ messages*age, data = social)
predicted <-
predict(mod1, newdata = data.frame(age      = rep(c(22:75), 2),
messages = rep(c("Control", "Neighbors"), each = 54)))
coef(mod1)
social <- read.csv("c:/users/jod.dbp/dropbox/undervisning/cbs/pmsd/lectures/prediction/social.csv")
social <- read.csv("~/dropbox/undervisning/cbs/pmsd/lectures/prediction/social.csv")
social$age <- 2008 - social$yearofbirth
social <- subset(social, messages == "Control" | messages == "Neighbors")
mod1 <- lm(100*primary2008 ~ messages + age, data = social)
mod1
predicted <-
predict(mod1, newdata = data.frame(age      = rep(c(22:75), 2),
messages = rep(c("Control", "Neighbors"), each = 54)))
plot(22:75, predicted[1:54], type = "n", ylim = c(10, 50), ylab = "Turnout", xlab = "Age", main = "Neighbor treatment over age")
lines(22:75, predicted[1:54], col = "Blue", lwd = 2)
lines(22:75, predicted[55:108], col = "Red", lwd = 2)
legend(x = 22, y = 50, legend = c("Neighbors", "Control"), lty = c(1,1), col = c("Red", "Blue"), lwd = 2, bty ="n")
lines(c(50,50), predicted[c(83, 93)])
lines(c(50,60), predicted[c(93, 93)], lty = 2)
text(y = predicted[93] - 2, x = 45, paste(10*round(coef(mod1)[3], 3), "(10*age) ="))
lines(c(60,60), predicted[c(39, 29)])
lines(c(50,60), predicted[c(29, 29)], lty = 2)
text(y = predicted[39] - 2, x = 65, paste("=", 10*round(coef(mod1)[3], 3), "(10*age)"))
lines(c(30,30), predicted[c(9, 63)])
text(x = 36, y = mean(predicted[c(9, 63)]) + 3, paste("=", round(coef(mod1)[2], 2), "(Neighbors)"))
summary(lm(100*primary2008 ~ age*messages, data = social))
mod1 <- lm(100*primary2008 ~ messages*age, data = social)
predicted <-
predict(mod1, newdata = data.frame(age      = rep(c(22:75), 2),
messages = rep(c("Control", "Neighbors"), each = 54)))
plot(22:75, predicted[1:54], type = "n", ylim = c(10, 50), ylab = "Turnout", xlab = "Age", main = "Neighbor treatment over age")
lines(22:75, predicted[1:54], col = "Blue", lwd = 2)
lines(22:75, predicted[55:108], col = "Red", lwd = 2)
legend(x = 22, y = 50, legend = c("Neighbors", "Control"), lty = c(1,1), col = c("Red", "Blue"), lwd = 2, bty ="n")
lines(c(25,25), predicted[c(4, 58)])
lines(c(70,70), predicted[c(49, 103)])
text(x = 30, y = mean(predicted[c(4, 58)]) + 3,
paste("=", round(coef(mod1)[2] + 25*coef(mod1)[4], 2), "(Neighbors)"))
text(x = 35, y = mean(predicted[c(4, 58)]) + 3,
paste("=", round(coef(mod1)[2] + 25*coef(mod1)[4], 2), "(b2 + 25*b4)"))
text(x = 35, y = predicted[4],
paste("=", round(coef(mod1)[2] + 25*coef(mod1)[4], 2), "(b2 + 25*b4)"))
text(x = 35, y = predicted[4],
paste("=", round(coef(mod1)[2] + 25*coef(mod1)[4], 2), "(b2 + 25*b4)"))
text(x = 45, y = predicted[103],
paste(round(coef(mod1)[2] + 25*coef(mod1)[4], 2), "(b2 + 70*b4) = "))
text(x = 65, y = predicted[103],
paste(round(coef(mod1)[2] + 25*coef(mod1)[4], 2), "(b2 + 70*b4) = "))
text(x = 60, y = predicted[103],
paste(round(coef(mod1)[2] + 25*coef(mod1)[4], 2), "(b2 + 70*b4) = "))
plot(22:75, predicted[1:54], type = "n", ylim = c(10, 50), ylab = "Turnout", xlab = "Age", main = "Neighbor treatment over age")
lines(22:75, predicted[1:54], col = "Blue", lwd = 2)
lines(22:75, predicted[55:108], col = "Red", lwd = 2)
legend(x = 22, y = 50, legend = c("Neighbors", "Control"), lty = c(1,1), col = c("Red", "Blue"), lwd = 2, bty ="n")
lines(c(25,25), predicted[c(4, 58)])
lines(c(70,70), predicted[c(49, 103)])
text(x = 35, y = predicted[4],
paste("=", round(coef(mod1)[2] + 25*coef(mod1)[4], 2), "(b2 + 25*b4)"))
text(x = 60, y = predicted[103],
paste(round(coef(mod1)[2] + 25*coef(mod1)[4], 2), "(b2 + 70*b4) = "))
marginal <- predicted[55:108] - predicted[1:54]
plot(22:75, marginal[1:54], type = "n", ylim = c(0, 20),
ylab = "Marginal effect \non turnout", xlab = "Age", main = "Neighbor treatment over age")
plot(22:75, marginal[1:54], type = "n", ylim = c(0, 20),
ylab = "Marginal effect on turnout", xlab = "Age", main = "Neighbor treatment over age")
lines(22:75, marginal[1:54], col = "Blue", lwd = 2)
par(mfrow = c (1,2))
plot(22:75, predicted[1:54], type = "n", ylim = c(10, 50), ylab = "Turnout", xlab = "Age", main = "Neighbor treatment over age")
lines(22:75, predicted[1:54], col = "Blue", lwd = 2)
lines(22:75, predicted[55:108], col = "Red", lwd = 2)
legend(x = 22, y = 50, legend = c("Neighbors", "Control"), lty = c(1,1), col = c("Red", "Blue"), lwd = 2, bty ="n")
lines(c(25,25), predicted[c(4, 58)])
lines(c(70,70), predicted[c(49, 103)])
text(x = 35, y = predicted[4],
paste("=", round(coef(mod1)[2] + 25*coef(mod1)[4], 2), "(b2 + 25*b4)"))
text(x = 60, y = predicted[103],
paste(round(coef(mod1)[2] + 25*coef(mod1)[4], 2), "(b2 + 70*b4) = "))
marginal <- predicted[55:108] - predicted[1:54]
plot(22:75, marginal[1:54], type = "n", ylim = c(0, 20),
ylab = "Marginal effect on turnout", xlab = "Age", main = "Neighbor treatment over age")
lines(22:75, marginal[1:54], col = "Blue", lwd = 2)
mod1 <- lm(primary2008 ~ messages*age, data = social)
mod3 <- lm(primary2008 ~ age + I(age^2), data = social)
summary(mod3)
predicted <- predict(mod3, newdata = data.frame(age = 22:75))
predicted
mod3 <- lm(primary2008 ~ age + I(age^2), data = social)
predicted <- predict(mod3, newdata = data.frame(age = 22:75))
predicted
mod3 <- lm(primary2008 ~ age + I(age^2), data = social)
mod3
predicted <- predict(mod3, newdata = data.frame(age = 22:75))
plot(22:75, predicted)
predicted <- predict(mod3, newdata = data.frame(age = 22:85))
plot(22:85, predicted)
plot(22:90, predicted, type = "n", ylim = c(10, 50), ylab = "Turnout", xlab = "Age", main = "Correlation between turnout and age")
plot(22:90, predicted, type = "n", ylim = c(10, 50), ylab = "Turnout", xlab = "Age", main = "Correlation between turnout and age")
predicted <- predict(mod3, newdata = data.frame(age = 22:90))
plot(22:90, predicted, type = "n", ylim = c(10, 50), ylab = "Turnout", xlab = "Age", main = "Correlation between turnout and age")
lines(22:90, predicted)
lines(22:90, predicted, col = "blue")
plot(22:90, predicted, type = "n", ylim = c(10, 50), ylab = "Turnout", xlab = "Age", main = "Correlation between turnout and age")
lines(22:90, predicted, col = "blue")
lines(22:90, predicted, col = "blue", lwd = 2)
plot(22:90, predicted, type = "n", ylim = c(10, 50), ylab = "Turnout", xlab = "Age", main = "Correlation between turnout and age")
lines(22:90, predicted, col = "blue", lwd = 2)
lines(x = 22:90, y = 100*predicted, col = "blue", lwd = 2)
x <- rnorm(n = 20, mean = 0, sd = 1)
y <- rnorm(n = 20, mean = 0, sd = 1)
z <- rnorm(n = 20, mean = 0, sd = 1)
summary(lm(y~x))$r.squared
set.seed(12345)
x <- rnorm(n = 15, mean = 0, sd = 1)
y <- rnorm(n = 15, mean = 0, sd = 1)
z <- rnorm(n = 15, mean = 0, sd = 1)
summary(lm(y~x))$r.squared
set.seed(12345)
x <- rnorm(n = 20, mean = 0, sd = 1)
y <- rnorm(n = 20, mean = 0, sd = 1)
z <- rnorm(n = 20, mean = 0, sd = 1)
summary(lm(y~x))$r.squared
summary(lm(y ~ x + z))$r.squared
str(mod3)
str(summary(mod3))
library(car)
install.packages("car")
library(car)
names(social)
table(social$hhsize)
scatter3d(primary2008 ~ age + hhsize, data = social)
library(rgl)
install.packages("rgl")
library(rgl)
library(rgl)
scatter3d(primary2008 ~ age + hhsize, data = social)
install.packages("plot3D")
library(plot3D)
scatter3d(x = social$age, y = social$primary2008, z =  social$hhsize)
detach(car)
library(plot3D)
social <- read.csv("c:/users/jod.dbp/dropbox/undervisning/cbs/pmsd/lectures/prediction/social.csv")
social <- read.csv("~/dropbox/undervisning/cbs/pmsd/lectures/prediction/social.csv")
social$age <- 2008 - social$yearofbirth
social <- subset(social, messages == "Control" | messages == "Neighbors")
mod1 <- lm(100*primary2008 ~ messages + age, data = social)
mod1
predicted <-
predict(mod1, newdata = data.frame(age      = rep(c(22:75), 2),
messages = rep(c("Control", "Neighbors"), each = 54)))
plot(22:75, predicted[1:54], type = "n", ylim = c(10, 50), ylab = "Turnout", xlab = "Age", main = "Neighbor treatment over age")
lines(22:75, predicted[1:54], col = "Blue", lwd = 2)
lines(22:75, predicted[55:108], col = "Red", lwd = 2)
legend(x = 22, y = 50, legend = c("Neighbors", "Control"), lty = c(1,1), col = c("Red", "Blue"), lwd = 2, bty ="n")
lines(c(50,50), predicted[c(83, 93)])
lines(c(50,60), predicted[c(93, 93)], lty = 2)
text(y = predicted[93] - 2, x = 45, paste(10*round(coef(mod1)[3], 3), "(10*age) ="))
lines(c(60,60), predicted[c(39, 29)])
lines(c(50,60), predicted[c(29, 29)], lty = 2)
text(y = predicted[39] - 2, x = 65, paste("=", 10*round(coef(mod1)[3], 3), "(10*age)"))
lines(c(30,30), predicted[c(9, 63)])
text(x = 36, y = mean(predicted[c(9, 63)]) + 3, paste("=", round(coef(mod1)[2], 2), "(Neighbors)"))
summary(lm(100*primary2008 ~ age*messages, data = social))
mod1 <- lm(100*primary2008 ~ messages*age, data = social)
predicted <-
predict(mod1, newdata = data.frame(age      = rep(c(22:75), 2),
messages = rep(c("Control", "Neighbors"), each = 54)))
par(mfrow = c (1,2))
plot(22:75, predicted[1:54], type = "n", ylim = c(10, 50), ylab = "Turnout", xlab = "Age", main = "Neighbor treatment over age")
lines(22:75, predicted[1:54], col = "Blue", lwd = 2)
lines(22:75, predicted[55:108], col = "Red", lwd = 2)
legend(x = 22, y = 50, legend = c("Neighbors", "Control"), lty = c(1,1), col = c("Red", "Blue"), lwd = 2, bty ="n")
lines(c(25,25), predicted[c(4, 58)])
lines(c(70,70), predicted[c(49, 103)])
text(x = 35, y = predicted[4],
paste("=", round(coef(mod1)[2] + 25*coef(mod1)[4], 2), "(b2 + 25*b4)"))
text(x = 60, y = predicted[103],
paste(round(coef(mod1)[2] + 25*coef(mod1)[4], 2), "(b2 + 70*b4) = "))
marginal <- predicted[55:108] - predicted[1:54]
plot(22:75, marginal[1:54], type = "n", ylim = c(0, 20),
ylab = "Marginal effect on turnout", xlab = "Age", main = "Neighbor treatment over age")
lines(22:75, marginal[1:54], col = "Blue", lwd = 2)
mod3 <- lm(primary2008 ~ age + I(age^2), data = social)
predicted <- predict(mod3, newdata = data.frame(age = 22:90))
plot(22:90, predicted, type = "n", ylim = c(10, 50), ylab = "Turnout", xlab = "Age", main = "Correlation between turnout and age")
lines(x = 22:90, y = predicted, col = "blue", lwd = 2)
library(plot3D)
scatter3D(x = social$primary2008, y = social$hhsize, z = social$age)
graphics.off()
data <- read.dta("turnout0513.dta") #data indl?ses
library(foreign)
graphics.off()
data <- read.dta("turnout0513.dta") #data indl?ses
data <- read.dta("~/dropbox/undervisning/cbs/pmsd/lectures/prediction/turnout0513.dta") #data indl?ses
data <- read.dta("~/dropbox/undervisning/cbs/pmsd/lectures/prediction/turnout0513.dta") #data indl?ses
names(data)
scatter3D(x = data$turnout05, y = data$turnout05, z = data$turnout05)
pt(1.64, df = 14)
1 - 2*pt(1.64, df = 14)
2*(1 - *pt(1.64, df = 14))
2*(1 - pt(1.64, df = 14))
2*(1 - pt(1.83, df = 14))
2*(1 - pt(1.27, df = 14))
2*(1 - pt(1.33, df = 14))
qnorm(0.975)
pnorm(qnorm(0.975))
STAR <- read.csv("~/dropbox/undervisning/cbs/pmsd/lab_session/data/STAR.csv")
7778
4500000*.1/136
4500000*.1/360
4500000*.05/360
4500000*.0295/360
9498/ 2 - 3000
9498/ 2 - 3000 - (1732.51 + 589 + 461) / 2
9498/ 2 - 3000 - (1732.51 + 589 + 466.66 + 392.7) / 2
11/21
12/22
3600*20
60000*0.02
4508 - 3861
4608*0.8
4608*0.2 - 100
4608*0.2
4608*0.2 - 100
1000 - 4608*0.2 - 100
4608*0.2 - 100
4608*0.2 - 100 - 1000
1000 - 4608*0.2 - 100
1000 - 4608*0.2 + 100
setwd(~/dropbox/own projects/turnout instrument)
setwd("~/dropbox/own projects/turnout instrument")
setwd("~/dropbox/own projects/turnout instrument/dataan/replication_files")
setwd("~/dropbox/own projects/turnout instrument/dataan/rd_all/replication_files")
# This files assess the robustness of Dahlgaard (2018)
# by varying the bandwidth from three to 183 days
load("agg_day_year_all.rdata")
load("n_day_year_all.rdata")
## commented out for replication files
#  install.packages("xtable")
#  install.packages("ggplot2")
#  install.packages("rmeta")
#  install.packages("reshape2")
#  install.packages("dplyr")
#  install.packages("rdrobust")
library(xtable)
library(ggplot2)
library(dplyr)
library(rmeta)
library(rdrobust)
library(reshape2)
## Define treatment dummy and assign weights
data_day_year <-
data_day_year %>%
ungroup() %>%
mutate(treated = days > 0)
#create large data set
data <- left_join(data_day_year, data_n_year,
by = c("days", "year"))
data_close <-
data %>%
mutate(voted     = round(par_vote*n),
abstained = round((1-par_vote)*n),
treated   = days > 0) %>%
select(year, days, treated, voted, abstained)
data_voted <-
as.data.frame(lapply(data_close,
function(x,p) rep(x,p),
data_close[["voted"]])) %>%
mutate(voted = 1)
data_abstained <-
as.data.frame(lapply(data_close,
function(x,p) rep(x,p),
data_close[["abstained"]])) %>%
mutate(voted = 0)
data_master <-
rbind(data_voted, data_abstained)
## Loop over bandwidths
effect <- data_frame(yr2009 = rep(NA, 183),
yr2013 = rep(NA, 183),
yr2014 = rep(NA, 183),
yr2015 = rep(NA, 183),
pooled = rep(NA, 183),
place  = 1:183,
rd     = "Conventional RD")
se <- data_frame(yr2009 = rep(NA, 183),
yr2013 = rep(NA, 183),
yr2014 = rep(NA, 183),
yr2015 = rep(NA, 183),
pooled = rep(NA, 183),
place  = 1:183,
rd     = "Conventional RD")
effect_robust <- data_frame(yr2009 = rep(NA, 183),
yr2013 = rep(NA, 183),
yr2014 = rep(NA, 183),
yr2015 = rep(NA, 183),
pooled = rep(NA, 183),
place  = 1:183,
rd     = "Robust RD")
se_robust <- data_frame(yr2009 = rep(NA, 183),
yr2013 = rep(NA, 183),
yr2014 = rep(NA, 183),
yr2015 = rep(NA, 183),
pooled = rep(NA, 183),
place  = 1:183,
rd     = "Robust RD")
years <- c(2009, 2013, 2014, 2015)
for (i in 3:183){
for(j in 1:4){
sub_data <-
data_master %>%
filter(year == years[j])
rd                  <- rdrobust(y = sub_data$voted, x = sub_data$days, h = i)
effect[i, j]        <- rd$coef[1]
se[i, j]            <- rd$se[1]
effect_robust[i, j] <- rd$coef[3]
se_robust[i, j]     <- rd$se[3]
}
effect[i, 5] <- meta.summaries(as.numeric(effect[i, 1:4]), as.numeric(se[i, 1:4]))$summary
se[i, 5]     <- meta.summaries(as.numeric(effect[i, 1:4]), as.numeric(se[i, 1:4]))$se.summary
effect_robust[i, 5] <-
meta.summaries(as.numeric(effect_robust[i, 1:4]), as.numeric(se_robust[i, 1:4]))$summary
se_robust[i, 5]     <-
meta.summaries(as.numeric(effect_robust[i, 1:4]), as.numeric(se_robust[i, 1:4]))$se.summary
if(i == 10*floor(i/10))
print(i)
}
effect_joint <-
rbind(effect, effect_robust) %>%
melt(., id = c("place", "rd"))
se_joint <-
rbind(se, se_robust) %>%
melt(., id = c("place", "rd")) %>%
rename(se = value)
bw_data <-
data_frame(variable = c("yr2009",
"yr2013",
"yr2014",
"yr2015",
"pooled"),
bw = round(c(with(data_master %>% filter(year == 2009),
rdbwselect(y = voted, x = days))$bws[1],
with(data_master %>% filter(year == 2013),
rdbwselect(y = voted, x = days))$bws[1],
with(data_master %>% filter(year == 2014),
rdbwselect(y = voted, x = days))$bws[1],
with(data_master %>% filter(year == 2015),
rdbwselect(y = voted, x = days))$bws[1],
NA)))
plot_data <-
left_join(effect_joint, se_joint,
by = c("place", "rd", "variable")) %>%
mutate(low = value + qnorm(0.025)*se,
hi  = value + qnorm(0.975)*se) %>%
left_join(., bw_data, by = "variable")
labels <- c(yr2009 = "2009",
yr2013 = "2013",
yr2014 = "2014",
yr2015 = "2015",
pooled = "Pooled")
plot <-
ggplot(data = plot_data,
aes(x = place, y = value)) +
geom_line() +
facet_grid(variable ~ rd,
labeller = labeller(variable = labels)) +
geom_line(aes(x = place, y = low), alpha = 0.5, linetype = 2) +
geom_line(aes(x = place, y = hi), alpha = 0.5, linetype = 2) +
theme_classic() +
xlab("Bandwidth in days") +
theme(axis.title = element_text(size = 14)) +
geom_vline(aes(xintercept = bw), alpha = 0.5, linetype = 2) +
geom_hline(yintercept = 0, alpha = 0.3, linetype = 3) +
scale_y_continuous("Effect in percentage points", limits = c(-0.15,0.15), breaks = seq(-0.1, 0.1, by = 0.05)) +
theme(axis.text = element_text(size = 14),
axis.title = element_text(size = 14),
legend.text = element_text(size = 14),
legend.title = element_text(size = 14))
plot
